﻿
class caja
{
private:

public:
	double lbox[3];//Lado de la caja
	particula *part; //Array de particulas
	int npart;//Numero de particulas
	double phi;//Concentracion en volumen
	double rho;//Concentracion en masa
	double temp;//Temperatura de caja
	double vol;//Volumen
	double rcut;//Distancia de corte de interaccion
	///Para saber tanto por ciento de movimientos aceptados
	double ntry;//Numero de intentos
	double nok; //Movimientos aceptados
	double u0;
	double dr[3];
	///Parametros lara la link-cell (vecinos)
	int ncell;//Numero de celdas 
	int ncell2;//ncell^2
	int ncell3;//ncell^3
	double lgrid;//Lado de celda
	/*head[i] contiene el primer indice de particula contenido en celda i
	head[i]=-1 si no hay particulas en celda i*/
	int *head;
	//Contiene la lista de particulas enlazadas
	int *link;
	/*Ejemplo:
	celda i
	p1=head[i]
		- p1 es el indice de la primera particula en celda i. 
		- si no hay particulas en celda i, p1=-1
	p2 = link[p1]
		- p2 es el indice de la segunda particula en celda i
		- p2=-1 si no hay mas particulas en la celda
	etc...
	Por lo tanto, en la celda i habrá
	[p1,p2,..., pk] particulas	
	*/
	//Energias
	double epoten;//Energia potencia por particula en la caja
	double virial;
	double ekin;//Energia cinetica por particula en la caja
	//Control de la velocidad (termostato)
	int nacc;///Paso de la simulacion
	double delta;//Paso de longitud
	parametros param;
	double u30;
	double temp0;

	caja()
	{

	}

	caja(parametros par)
	{
		///Empleamos paso de tiempo Dt para generar las posiciones posOld
		///de las particulas empleando las velocidades
		param=par;
		nok=0;
		ntry=0;
		if(param.rest==0)
		{
			//Iniciamos la simulacion desde el principio
			temp=1.0;
			temp0=1.0;
			u0=param.u0;
			u30=3.0*u0;
			npart=param.npart;
			temp0 =1.0;
			phi=param.phi;
			rcut = param.rcut;
			rho = phi*6.0/pi;
			delta=param.Dr;
			double bx=pow((4./3.)*pi*npart/phi,1./3.);

			for(int i=0;i<3;i++)
				lbox[i]= bx;
			part = new particula[npart];
			link = new int[npart];
			vol=lbox[0]*lbox[1]*lbox[2];

			ncell=int(lbox[0]/rcut);
			ncell2=ncell*ncell;
			ncell3=ncell*ncell*ncell;
			lgrid= lbox[0]/double(ncell);
			cout<<"Numero de celdas: "<<ncell<<endl;
			cout<<"Lado de celda: "<<lgrid<<endl;
			head = new int[ncell*ncell*ncell];

			SetMeshPart();
		}
		else
		{
			//Reiniciamos una simulacion desde archivo existente
			loadConfig();
		}
	}

	caja(int np, double xT, double xphi, double xrcut, double del)
	{
		///Empleamos paso de tiempo Dt para generar las posiciones posOld
		///de las particulas empleando las velocidades
		ntry=0;
		nok=0;
		delta=del;
		npart=np;
		temp =xT;
		phi=xphi;
		rcut = xrcut;
		rho = phi*6.0/pi;
		double bx = pow(3.1416*npart/(6.*phi),1./3.);

		for(int i=0;i<3;i++)
			lbox[i]= bx;
		part = new particula[npart];
		link = new int[npart];
		SetMeshPart();
		vol=lbox[0]*lbox[1]*lbox[2];

		ncell=int(lbox[0]/rcut);
		ncell2=ncell*ncell;
		ncell3=ncell*ncell*ncell;
		lgrid= lbox[0]/double(ncell);
		cout<<"Numero de celdas: "<<ncell<<endl;
		cout<<"Lado de celda: "<<lgrid<<endl;
		head = new int[ncell*ncell*ncell];
	}

	///Reinicia la simulacion desde un archivo
	caja(char filename[50], double del)
	{
		ntry=0;
		nok=0;
		delta=del;
		///Empleamos paso de tiempo Dt para generar las posiciones posOld
		///de las particulas empleando las velocidades
		ifstream MiArchivo (filename,ios::in);
		if (MiArchivo.is_open())
		{
			///Cargamos las caracteristicas de la caja
			MiArchivo>>npart>>phi>>temp>>rcut;
			rho = phi*6.0/3.1416;
			double bx = pow(3.1416*npart/(6.*phi),1./3.);
			for(int i=0;i<3;i++)
				lbox[i]= bx;
			part = new particula[npart];
			link = new int[npart];
			SetMeshPart();
			vol=lbox[0]*lbox[1]*lbox[2];

			ncell=int(lbox[0]/rcut);
			ncell2=ncell*ncell;
			ncell3=ncell*ncell*ncell;
			lgrid= lbox[0]/double(ncell);
			cout<<"Numero de celdas: "<<ncell<<endl;
			cout<<"Lado de celda: "<<lgrid<<endl;
			head = new int[ncell*ncell*ncell];
			///Cargamos las posiciones de las particulas
			for(int i=0;i<npart;i++)
			{
				for(int j=0;j<3;j++)
				{
					MiArchivo>>part[i].pos[j];
					//Posiciones anteriores de las particulas
					part[i].posOld[j]=part[i].pos[j];
				}
			}
			MiArchivo.close();
		}

	}

	~caja()
	{
		delete [] part;
	}

	///Imprime la configuracion
	void printConfig(char namefile[30]);
	///Coloca las particulas de forma aleatoria en una caja
	void SetRandomPart(void);
	///Coloca las particulas en red
	void SetMeshPart(void);
	///Calcula la linkCell
	void linkCell(void);
	///Aplica condiciones contorno periodicas a particula ind
	void boundary(int ind);
	///Calcula la distancia cuadratica entre particulas
	double distancia2(int ind1,int ind2);
	///Calcula la distancia entre particulas
	double distancia(int ind1,int ind2);
	///Imprime el archivo XYZ de la configuracion de particulas
	void printXYZfile(char namefile[30]);
	//Selecciona una particula al azar y la desplaza
	void randomMovement(void);
	//Inicializa la energia de cada particula
	void calcEnergia(void);
	//Calcula la energia de la particula ind
	double calcEnergiaParticle(int ind);
	///Elimina a una particula de una celda
	void removeParticleFromCell(int ind);
	///Inserta particula en una celda
	void insertParticleInCell(int icell, int ind);
	void loadConfig(void);
	double etotal(void);
	double caja::calcEnergiaParticleLCL(int index);

};

double caja::calcEnergiaParticle(int index)
{
	int c, c1;
	int iii,jjj;
	int shift[3];
	double dist, dist2, dist3, potdm, ePoten;
	double rcut2=rcut*rcut;
	double f2, Co, Co2;
	int count=0;

	ePoten=0.;

	for(int iii=0;iii<npart-1;iii++)
	{
		for(int jjj=iii+1;jjj<npart;jjj++)
		{
			///Calculamos la distancia con condiciones periodicas entre iii y jjj
			dist2= distancia2(iii,jjj);
			//Si la distancia es menor que el cutoff calculamos la interaccion
			dist=sqrt(dist2);
			dist3=dist*dist*dist;
			count++;
			if (dist<4.0)
			{
				//Para el potencial
				f2=0.00375*u0*exp(-50.*(dist-2.));
			}
			else
			{
				f2=0.0;
			}

			///Calculamos la fuerza entre iii y jjj
			Co=dr[2]/dist;
			Co2=Co*Co;																	

			potdm = f2+(u0)*(1.0-3*Co2)/(dist3);
			//potdm=0.0;
			//Acualizamos la energia potencial
			ePoten+=potdm;
		}
	}

	return ePoten;
}

double caja::etotal(void)
{
	double etotal=0.0;
	double emax=-1.0e8;
	double emin=1.0e8;
	for(int i=0;i<npart;i++)
	{
		if(part[i].energia>emax)
			emax=part[i].energia;
		if(part[i].energia<emin)
			emin=part[i].energia;
		etotal+=part[i].energia;
	}
	epoten=etotal/double(npart);
	return etotal;
}

void caja::randomMovement(void)
{
	double alea = (double(rand())/RAND_MAX);
	int index = int(alea*npart);
	if (index==npart) index--;
	double posR;
	int mc[3], mcOld[3];
	int icellOld;
	int icellNew;
	double energia;

	//Actualizamos el numero de intentos
	ntry++;
	//Celda en la que se encontraba la particula
	icellOld=part[index].celda;
	//Desplazamos aleatoriamente la particula
	double enold = calcEnergiaParticle(index);
	for (int j=0;j<3;j++)
	{
		posR=(part[index].pos[j])+delta*(double(rand())/RAND_MAX-0.5);
		part[index].posOld[j]=part[index].pos[j];
		part[index].pos[j]=posR;
	}
	//Condiciones de contorno periodicas
	boundary(index);
	//Localizamos la celda en la que se encuentra la particula
	for(int j=0;j<3;j++)
	{
		mc[j] = int((part[index].pos[j]+0.5*lbox[0])/lgrid);
	}
	icellNew = int(mc[0]*ncell2+mc[1]*ncell+mc[2]);
	if(index==32 && icellNew==57)
	{
		int jjx=123;
	}
	
	if(icellNew!=icellOld)
	{
		//Actualizamos las coordenadas de la celda donde se encuentra la particula
		for(int j=0;j<3;j++)
		{
			mcOld[j]=part[index].celdaxyz[j];
			part[index].celdaxyz[j]=mc[j];
		}
		//La particula ha cambiado de celda
		//Sacamos a la particula de la celda en la que se encontraba
		removeParticleFromCell(index);
		//Introducimos a la particula en la nueva celda
		insertParticleInCell(icellNew, index);
		//Acualizamos el objeto 'particula' con la nueva celda
		part[index].celda=icellNew;		
	}

	///Calculamos la nueva energia de la particula
	energia = calcEnergiaParticle(index);

	//if(energia<=part[index].energia)
	if(energia<=enold)
	{
		//Actualizamos el numero de intentos satisfactorios
		nok++;
		//cout<<"aceptado"<<endl;
		//Aceptamos el movimiento y actualizamos la energia
		//part[index].energiaOld=part[index].energia;
		//part[index].energia=energia;
	}
	else
	{
		alea = (double(rand())/RAND_MAX);
		//double xx= -(energia-part[index].energia)/temp;
		double xx= -(energia-enold)/temp;
		double tt = exp(xx);
		if(alea<tt)
		{
			//Actualizamos el numero de intentos satisfactorios
			nok++;
			//cout<<"aceptado"<<endl;
			//Aceptamos el movimiento y actualizamos la energia
			//part[index].energiaOld=part[index].energia;
			//part[index].energia=energia;
		}
		else
		{
			//cout<<"rechazado"<<endl;
			//Rechazamos el movimiento y colocamos a la particula en su posicion anterior
			for(int j=0;j<3;j++)
			{
				part[index].pos[j]=part[index].posOld[j];
			}
			if(icellNew!=icellOld)
			{
				for(int j=0;j<3;j++)
					part[index].celdaxyz[j]=mcOld[j];
				//Tenemos que colocar de nuevo a la particula en la celda que le corresponde
				removeParticleFromCell(index);
				//Introducimos a la particula en la nueva celda
				insertParticleInCell(icellOld, index);
				part[index].celda=icellOld;
			}

		}
	}
}

double caja::calcEnergiaParticleLCL(int index)
{
	int c, c1;
	int iii,jjj;
	int shift[3];
	double dist, dist2, dist3, potdm, ePoten;
	double rcut2=rcut*rcut;
	double f2, Co, Co2;
	int count=0;

	ePoten=0.;

	c = part[index].celda;
	///Recorremos las particulas de esta celda y las celdas vecinas
	for (int ixred=part[index].celdaxyz[0]-1;ixred<=part[index].celdaxyz[0]+1;ixred++)
	{
		/*shift: condiciones de contorno periodicas tenidas en cuenta para
		localizar las celdas vecinas*/
		if (ixred < 0)
		{
			shift[0] = int(ixred+ncell);
		}
		else
		{
			if (ixred>=ncell)
			{
				shift[0] = int(ixred-ncell); 
			}
			else
			{
				shift[0] = int(ixred);
			}
		}
		for (int iyred=part[index].celdaxyz[1]-1;iyred<=part[index].celdaxyz[1]+1;iyred++)
		{
			if (iyred < 0)
			{
				shift[1] = int(iyred+ncell);
			}
			else
			{
				if(iyred>=ncell)
				{
					shift[1] = int(iyred-ncell);
				}
				else
				{
					shift[1] = int(iyred);
				}
			}
			for (int izred=part[index].celdaxyz[2]-1;izred<=part[index].celdaxyz[2]+1;izred++)
			{
				if (izred < 0)
				{
					shift[2] = int(izred+ncell);
				}
				else
				{
					if(izred>=ncell)
					{
						shift[2] = int(izred-ncell);
					}
					else
					{
						shift[2] = int(izred) ;
					}
				}
				///Calculamos el indice c1 de la celda vecina
				c1 = int(((shift[0]+ncell)%ncell)*ncell2+((shift[1]+ncell)%ncell)*ncell+((shift[2]+ncell)%ncell));
				///El indice index ahora es iii
				iii = index;

				///Buscamos la primera particula jjj de la celda c1
				jjj =head[c1];

				while (jjj != -1) ///jjj=-1->No hay particulas en la celda
				{
					if(iii!=jjj)
					{
						///Calculamos la distancia con condiciones periodicas entre iii y jjj
						dist2= distancia2(iii,jjj);
						//Si la distancia es menor que el cutoff calculamos la interaccion
						if (dist2 < rcut2 )
						{
							dist=sqrt(dist2);
							dist3=dist*dist*dist;
							count++;
							if (dist<4.0)
							{
								//Para el potencial
								f2=0.00375*u0*exp(-50.*(dist-2.));
							}
							else
							{
								f2=0.0;
							}

							///Calculamos la fuerza entre iii y jjj
							Co=dr[2]/dist;
							Co2=Co*Co;																	

							potdm = f2+(u0)*(1.0-3*Co2)/(dist3);
							//potdm=0.0;
							//Acualizamos la energia potencial
							ePoten+=potdm;							
						}
					}
						//Pasamos a la siguiente particula de la celda c1
						jjj = link[jjj];
				}
			}
		}
	}
	return ePoten;
}

void caja::SetRandomPart(void)
{
	double pos[3];

	for(int i=0;i< npart;i++)
	{
		for (int j=0;j<3;j++)
		{
			pos[j]=lbox[j]*double(rand())/RAND_MAX;
		}
		part[i].setPosition(pos);
		boundary(i);
	}
}

void caja::removeParticleFromCell(int ind)
{
	int icell = part[ind].celda;
	int i = head[icell];
	int j, iaux;
	if(i==ind)
	{
		j=link[i];
		head[icell]=j;
	}
	else
	{
		while (i!=-1)
		{
			iaux = link[i];
			if (iaux==ind)
			{j=link[iaux];}
			else
			{j=iaux;}
			link[i]=j;
			i=j;
		}
	}


}

void caja::insertParticleInCell(int icell, int ind)
{
	int iold;
	int i = head[icell];
	if (i==-1)
	{
		head[icell]=ind;
		link[ind]=-1;
	}
	else
	{
		while (i!=-1)
		{
			iold=i;
			i = link[i];
			if (i==-1)
			{
				link[iold]=ind;
				link[ind]=-1;
				i=-1;
			}
		}
	}
}

//Listas 

void caja::boundary(int ind)
{
	for(int i=0;i<3;i++)
	{
		part[ind].pos[i]-=lbox[i]*anint(part[ind].pos[i]/lbox[i]);
	}
}

double caja::distancia2(int ind1,int ind2)
{
	//Calcula la distancia entre las particulas
	double dis=0;
	for (int i=0; i<3;i++)
	{
		dr[i]=part[ind2].pos[i]-part[ind1].pos[i];
		dr[i]-=lbox[i]*anint(dr[i]/(lbox[i])); //Entero + proximo manteniendo signo 
		dis+=dr[i]*dr[i];
	}
	return	 dis;
}

double caja::distancia(int ind1,int ind2)
{
	return	 sqrt(distancia2( ind1,ind2));
}

void caja::printXYZfile(char namefile[30])
{
	char name[50];
	sprintf(name, "%s%s", namefile,".xyz");
	ofstream MiArchivo (name,ios::out|ios::app);
	MiArchivo.precision(5);
	if (MiArchivo.is_open())
	{
		MiArchivo << npart<<endl;
		MiArchivo << lbox[0]<<"  "<<nacc<<endl;
		for(int i=0;i<npart;i++)
		{
			MiArchivo<<setw(5)<<"O "<<setw(8)<<part[i].pos[0]<<"  "<<setw(8)<<part[i].pos[1]<<"  "<<setw(8)<<part[i].pos[2]<<endl;
		}
		MiArchivo.close();
	}
}

void caja::linkCell(void)
{
	int mc[3], c;
	//Reset the headers, head
	for (int i=0;i<ncell3;i++)
	{
		head[i] = -1;
	}

	///Scan atoms to construct headers, head, & linked lists, lscl
	for (int i=0;i<npart;i++)
	{
		///Vector cell index to which this atom belongs
		for (int j=0;j<3;j++)
		{
			mc[j] = int((part[i].pos[j]+0.5*lbox[0])/lgrid);
		}
		///Translate the vector cell index, mc, to a scalar cell index"""
		//
		c = int(mc[0]*ncell2+mc[1]*ncell+mc[2]);
		part[i].celda=c;
		part[i].celdaxyz[0]=mc[0];
		part[i].celdaxyz[1]=mc[1];
		part[i].celdaxyz[2]=mc[2];

		///Link to the previous occupant (or EMPTY if you're the 1st)
		link[i] = head[c];
		///The last one goes to the header"""
		head[c] = i;
	}

	int ntotal=0;
	for(int i=0;i<ncell3;i++)
	{
		int cual = head[i];
		int n=0;
		while(cual!=-1)
		{
			cual=link[cual];
			n++;
		}
		ntotal+=n;
	}
	cout<<"Numero total de particulas colocadas en Mesh "<<ntotal<<endl;
}

void caja::printConfig(char namefile[30])
{
	char name[50];
	sprintf(name, "%s%s%d%s", namefile,"_",nacc,".txt");
	ofstream MiArchivo (name,ios::out);
	MiArchivo.precision(4);

	

	if (MiArchivo.is_open())
	{
		MiArchivo.setf(ios::fixed);
		MiArchivo<<"Parametros de la simulacion"<<endl;
		MiArchivo.fill('*');
		MiArchivo << setw(122) << ""<<endl;
		MiArchivo.fill(' ');
		MiArchivo<<setw(10)<<"Npart"<<setw(10)<<"phi"<<setw(10)<<"rcut"<<setw(10)<<"Dr"<<setw(10)<<"restart"<<setw(10)<<"bx"<<setw(10)<<"by"<<setw(10)<<"bz"<<setw(10)<<"epoten"<<endl;
		MiArchivo.fill('-');
		MiArchivo << setw(122) << ""<<endl;
		MiArchivo.fill(' ');
		MiArchivo<<setw(10)<<npart<<setw(10)<<phi<<setw(10)<<rcut<<setw(10)<<param.Dr<<setw(10)<<param.rest<<setw(10)<<lbox[0]<<setw(10)<<lbox[1]<<setw(10)<<lbox[2]<<setw(10)<<epoten<<endl;
	}
	MiArchivo<<endl;
	MiArchivo<<"Posiciones"<<endl;
	MiArchivo.fill('*');
	MiArchivo << setw(122) << ""<<endl;
	MiArchivo.fill(' ');
	MiArchivo<<setw(20)<<"rx"<<setw(20)<<"ry"<<setw(20)<<"rz"<<endl;
	MiArchivo.fill('-');
	MiArchivo << setw(122) << ""<<endl;
	MiArchivo.fill(' ');
	for(int i=0;i<npart;i++)
	{
		for(int j=0;j<3;j++)
			MiArchivo<<setw(20)<<part[i].pos[j];
		MiArchivo<<endl;
	}
}

void caja::loadConfig(void)
{
	ifstream MiArchivo (param.iniFile,ios::in);
	char line[150];
	char dummy[20];
	if(MiArchivo.fail())
	{
		cout<<"No se pudo encontrar el archivo de reinicio "<<param.iniFile<<endl;
		exit(1);
	}
	if (MiArchivo.is_open())
	{
		MiArchivo.getline(line,150);
		MiArchivo.getline(line,150);
		MiArchivo.getline(line,150);
		MiArchivo.getline(line,150);

		MiArchivo>>param.npart>>param.phi>>param.rcut>>param.Dr>>param.rest>>lbox[0]>>lbox[1]>>lbox[2]>>dummy;
		MiArchivo.getline(line,150);
		MiArchivo.getline(line,150);
		MiArchivo.getline(line,150);
		MiArchivo.getline(line,150);
		MiArchivo.getline(line,150);
		MiArchivo.getline(line,150);

		//Construimos las cosas
		u0=param.u0;
		u30=3.0*u0;
		npart=param.npart;
		temp0 =1.0;
		temp=1.0;
		phi=param.phi;
		rcut = param.rcut;
		rho = phi*6.0/pi;
		delta=param.Dr;
		part = new particula[npart];
		link = new int[npart];
		vol=lbox[0]*lbox[1]*lbox[2];
		ncell=int(lbox[0]/rcut);
		ncell2=ncell*ncell;
		ncell3=ncell*ncell*ncell;
		lgrid= lbox[0]/double(ncell);
		cout<<"Numero de celdas: "<<ncell<<endl;
		cout<<"Lado de celda: "<<lgrid<<endl;
		head = new int[ncell*ncell*ncell];
		//Fin de la construccion
		
		//Importamos las posiciones y velocidades
		for(int i=0;i<npart;i++)
		{
			for(int j=0;j<3;j++)
			{
				double rr;
				MiArchivo>>rr;
				part[i].pos[j]=rr;
			}
			boundary(i);
		}

		MiArchivo.close();
	}
}

void caja::SetMeshPart(void)
{
	int nmax3=anint(pow(npart,1./3.));
	int ipart=0;
	double lg=lbox[0]/(pow(double(npart),(1./3.)));
	double alea;

	for(int i=1;i<=nmax3;i++)
	{
		for(int j=1;j<=nmax3;j++)
		{
			for(int k=1;k<=nmax3;k++)
			{
				alea=(double(rand())/RAND_MAX-0.5);
				part[ipart].pos[0]=-lbox[0]/2.+((2.0*double(i)-1.0)/2.)*lg+lg*(1./4.)*alea;
				alea=(double(rand())/RAND_MAX-0.5);
				part[ipart].pos[1]=-lbox[1]/2.+((2.0*double(j)-1.0)/2.)*lg+lg*(1./4.)*alea;
				alea=(double(rand())/RAND_MAX-0.5);
				part[ipart].pos[2]=-lbox[2]/2.+((2.0*double(k)-1.0)/2.)*lg+lg*(1./4.)*alea;
				ipart++;
			}
		}
	}
}